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I. Overall Objectives 

The principal objective of this project is to apply the method of full 
simulation to reacting turbulent flows. Full simulation has proven of great 
value as a complement to experiments for the study of nonreacting turbulent 
flows. It provides insight into the physics of turbulent flows and their 
modeling. It is natural to try to extend these methods to the simulation of 
reacting turbulent flows. Because this is one of the first attempts at this 
type of simulation, a subsidiary goal of this work is to demonstrate the 
feasibility of using simulation to study turbulent reacting flows. In addi- 
tion, we intend to show that such simulations can be used to provide physical 
insight into the nature of turbulent combustion and to provide data that will 
help to construct models that can be used in engineering simulations of turbu- 
lent reacting flows. 

Reacting flows are considerably more difficult to simulate than non- 
reacting flows. The presence of chemical reaction demands that conservation 
equations for the various chemical species be solved along with the equations 
describing the fluid mechanics. Thus, the simulation of combusting flows 
includes all of the difficulties found in nonreacting turbulent flows as well 
as an equally large set associated with computing chemically reacting flows. 
However, experimental measurements are also very difficult to make in reacting 
turbulent flows, so the benefits produced by a successful simulation are 
greater; this is a primary motivation -for pursuing simulations of this type. 

We first had to choose the problems to be considered; this choice is im- 
portant in any research effort, but even more important in this one, due to 
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the need for this work to demonstrate the feasibility of the method and to 
provide the basis for farther simulations of a more applied nature. Since the 
grant supports two student research assistants, two problems had to be selec- 
ted. The choice was based on the following criteria: (1) relative simplicity 
of both geometry and chemistry; (2) this is essential in order to provide fea- 
sibility of simulation in a reasonable time on current computers; (3) further- 
more, to allow careful checking of the results, we need to consider flows for 
which sufficient experimental data are available; (4) finally, we want prob- 
lems which, to the extent possible, have relevance to applications. In addi- 
tion, the problems should be as different from each other as possible. After 
considerable searching, we selected the reacting axisymmetric jet and turbu- 
lent flame propagation. We shall describe the work on each of these problems' 
separately below. 

II. Reacting Axisymmetric Jet 

A. Physical Description 

Shear is fundamental to the mixing of fuel and oxidant in nearly all 
combustors in which the reactants are not premixed. The presence of strong 
shear ensures the creation of turbulence that enhances the mixing process; 
mixing must take place before chemical reaction is possible. In most combus- 
tors, the flow configuration is essentially that of an axisymmetric jet; the 
external air stream may be stationary or moving at a slower speed than the 
fuel or fuel-air jet. For this reason, we have chosen to study the axisymmet- 
ric jet-diffusion flame as one of the initial cases. This flow has been the 
subject of recent experiments at Stanford, and our knowledge of this flow will 
be very valuable in helping to simulate it. 

The nonreacting axisymmetric jet has not yet been treated by full simula- 
tion. However, a related flow, the spatially developing mixing layer, has 
recently been simulated by one of the members of our group. This flow shares 
many of the features of the axisymmetric jet but differs in important ways 
from the time-developing mixing layers that have been simulated by a number of 
other authors. Several important problems had to be resolved in order to make 
this simulation. Firstly, the conditions at the upstream boundary of the com- 
putational domain have to be specified in full detail; these are not available 
from the experiment but are important because they exert a strong influence 
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upon the flow. This problem was dealt with by constructing conditions repre- 
sentative of a transitional mixing layer and investigating their influence on 
the flow. Secondly, the condition at the outflow boundary is also difficult 
to specify. A condition that does not affect the solution upstream of the 
boundary was found. Thirdly and finally, a numerical method for computing the 
solution with these boundary conditions was constructed. 

These methods will be extended to treat the axisymmetric jet. This will 
require the development of new numerical methods for dealing with problems in 
cylindrical geometry. For the mixing layer, careful experiments have shown 
that the two- and three-dimensional flows have different rates of mixing. 
This is not yet so clear for the jet, but there is reason to believe that it 
is true, and we shall adopt it as a working hypothesis. For this reason, we 
intend to do a relatively large number of two-dimensional simulations, but it 
is important that some three-dimensional cases be done so that comparisons of 
the two cases can be made. 

These arguments provide the basis for the directions in which we intend 
to go. First, it is necessary to develop the numerical methods for dealing 
with axisymmetric flows. Then, simulations of completely axisymmetric jets 
(one with no azimuthal dependence of any of the variables) will be made. 
These will include flows with combustion, which will be treated by a simple 
model, and comparison cases without chemical reaction. All of the initial 
cases will be completely axisymmetric in order to reduce the computer— time 
requirements and turnaround time. Studies will be made of the effects of the 
heat-release parameter and the velocity ratio in the jet on the nature of the 
flame and its properties. When studies of the axisymmetric case are complete, 
a few carefully selected nonaxisymmetric cases will be run with a view to 
studying the differences between them and the completely axisymmetric cases. 

We shall begin with a description of the numerical method we have devel- 
oped. 
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B. Numerical Method for Axisymmetric Flows 


The Hankel transform is ideal for numerically solving problems in axisym- 
metric geometry and was actively considered. The method was found useful for 
problems in which the solution decays rapidly with radial distance. Several 
numerical versions of the Hankel transform were investigated. All of these 
are approximations to the exact integral transform, but, unfortunately, none 
of them is invertible, in the sense that one cannot retrieve the original data 
by a second transform. Furthermore, the convergence properties of Hankel 
series are less than desirable. 

We then decided to restrict our attention to finite geometry, in which 
the range of the radial variable is finite. This permits use of several fam- 
ilies of orthogonal polynomials in the radial variable. 

Some desirable properties of the polynomials are: 

• Orthogonality over (0,1) domain. Any domain, say, 0 _< r _<_ r Q , can be 
readily transformed to this domain. 

• Solutions of a Sturm-Liouville ordinary differential equation that is 
singular at r = 0 and r * 1. This enables application of arbitrary 
boundary conditions without sacrifice of spectral accuracy. For the 
axisymmetric problem, the radial derivative at r 31 0 should vanish, but 
this is unnecessary for the nonaxisymmetric case. 

• High-order members should have their zeroes clustered in regions where 
high resolution is desired. In our problem, this could be done at one 
station. However, at different downstream locations, the region requir- 
ing resolution shifts, and the polynomials may no longer be optimal. 

• Polynomials should admit a fast transform, enabling easy transformation 
in both directions. By fast, we mean methods requiring less than 0(n J ) 
operations. At the moment, such methods are not available for our choice 
of polynomials. 

The basis functions chosen are the normalized Jacobi polynomials. The 
standard Jacobi polynomials are defined for r e (0,1) and denoted by 
G n (p>q>r). G n is a polynomial of degree n with weight function 

w(r) - (l-r) P_<1 r ql ; p-q > 1, q>0 (1) 
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We choose p = q = 0.5 so that 


w(r) * l//r . (2) 

From here on, we shall call these particular polynomials G n (r). The differ- 
ential equation satisfied by them is given in [1]. 

r(l-r) G"(r) + (0.5 - 5r) G'(r) + n(n+0.5) G (r) - 0 (3) 

n n n 

or, in Sturm-Liouville form, 

[/7 (1-r) G' ( r ) ] ' + n(n+0. 5) G (r)//7 - 0 (4) 

n n 

Note that, since P(r) = /r (1-r) is zero at r ■ 0 and 1, the ODE (4) is 
singular at these points. 

We chose as the basis functions: 

H (r) = G (r)//a~ (5) 

n n n 

so that H n (r) satisfies the orthonormality relation, 

H (r) H (r) [1//7] dr = 6 (6) 

n m nm 

The ODE (4) is also satisfied by Hn(r). The recursion relation satisfied by 

the H is: 
n 

rH (r) = H . (r) /b(n+l) + H (r) a(n) + H . (r) /b(n) (7) 

n n+1 n n-i 

Using the expression for the derivative of P n ^ a ’^(r) [1], we can derive the 

following relationship satisfied by H^(r): 

r(l-r)H'(r) = H ,,(r) d(n) + H (r) e(n) + H .(r)f(n) (8) 

n n+1 n n— l 

where 



d(n) = - n/b(n+l) 


e(n) 


2 


n 

(2n - 0.5) 


na(n) 
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f( 


■> - {?£ 


(2n-l) (n-0.5) 


1) (4n-3) (2n-0.5) 


- nb( 


n)]^- 

/b(n) 


The zeros can be found by the method of Welsch and Golub [2]. 

We shall approximate a function f(r) by its values at N points Tq, 
r l» '•* r N-l’ we ca -*-l t * le approximant f^(r) an< * it by: 

N-l 

f N (r i ) ’ 2 H j <r l> • <9) 

j-0 

Requiring that f N (r^) “ f(r^), the exact value, leads to a system of N 

A A A 

equations for the N unknowns f Q , f^, ... This may be represented as 


A f 


In 


(l) 


where A is a full matrix of size N x N, whose elements are 


ij 


(r £ ) , 


The inverse of (10) is 


i - 0, 1, ... N-l, j - 0, 1 (N-l) 


f 



(ID 


The process of finding _f from f^ is called the forward Jacobi transform; 
the inverse is called the (backward) Jacobi transform. In general, these 
transforms require operation counts of 0(N ) as noted earlier. 


We have not yet specified the N collocation points. If these are cho- 
sen to be the zeroes of the N tll -order polynomial, certain attractive proper- 
ties result. 


If, in Eq. (9), we remove the subscript i from r on both sides, we 
get a continuous function of r, an interpolant of the data provided by 
f^(r^). This interpolant is at most a polynomial of degree (N-l) in r. 
We want to use it to estimate derivatives. 


Differentiating both sides of the interpolant with respect to 
using the derivative relationship (8), we have 
N-l 

df N /dr - ^ fj [d(j )Hj +1 (r) + e(j)Hj(r) + f(j)Hj_ 1 (r)]/r(l-r) « 

j-o 


r and 


( 12 ) 
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Since this is a polynomial of degree no greater than (N-2) in r, we may 
write 

N-l 

df jj/ dr = ^ Sj H j( r ) (13) 

j=0 

Clearly, gfj-i must be zero. Our aim is to relate the j* to the f_. 
Multiplying (13) by r(l-r) and using (12), 


r(l-r) df N dr 


N-l A i— 

f j |^ (j )H j+l (r) + + f <J> H j-l (r) 

N-l 

£ r(l-r) g.H (r) 
j=0 3 3 


(14) 


By repeated use of the recurrence relation (7) followed by index shifts, we 
get 

Z [f d(J-l) + fje(J) + f j+1 f(j+Dl Hj (r) + f N _ 1 d(N-l)H N (r) 

- E [g J _ 2 {-V r b(j-l) /b(j)> + g j _ 1 /b(j) {l-a(j) - a(j-l)} + 

+ g.(a(j) ~ b(j+l) - a(j) 2 - b( j )} + g ./b(j+l) {1 - a(j+l) - 

j . 3 (15) 

- a(j» + g j+2 (- /b(j+2) /b(j+l)}] H (r) + 

+ H N (r) {g^ /bOO (l-a(N) - a(N-l)) - g N _ 2 /b(N) /b(N-l)} + 

+ H N+l (r) { " %-l •'MN+T)} • 

A A 

It is assumed that, when the index on f or g falls outside the range 
(0,N-1), it is set to zero. Comparing coefficients of Jacobi polynomials of 
the same order, we get 

A 

Coefficient of H N+1 (r) -► g JJ _ 1 » 0 , 

Coefficient of H N (r) * g N _ 2 * f N _jd(N-l)/{-/b(N) /b(N-l)} 

• A 

Coefficient of H 2 (r) g^ can be obtained • 
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Thus the process of obtaining the g , the Jacobi coefficients of df/dr is 

J ~ 

quite efficient; df/dr can then be computed by inverting j* • 


Using similar ideas, it is possible to obtain other derivatives of f 
including d(rdf jj/dr)/dr, which occurs in the Navier-Stokes equations. 

Next, we shall consider the application of these methods to the solution 
of first- and second-order ODEs. 


Example 2 : 

Consider solving a first-order ODE, 


df (r) 
dr 


e(r) 


f(r) = f 


(16a) 

(16b) 


We shall use the spectral method with an (N+l) point approximation to f(r). 
The additional point enables us to satisfy the initial condition (16b). Thus, 
let us seek a solution of the form 

N 

f N-H (r) * 2 'jV r) <17> 

j-0 


Let us define the inner product of two functions: 

r 1 

(f i < 

‘ o 


dr 


(18) 


(N + l) equations for the (N+l) unknowns f . . • f^ are needed. N 


equations are obtained from the N inner products. 


(Lf N+1 ,H fc (r) r(l-r)) = (9(r), H k (r) r(l-r)) , 

k = 0, 1, 2, ... , (N-l ) . 

The last equation is derived from the initial condition: 

N 


09) 


N+l 


(0) = £ f H (0) = f 

j“0 J ( 


( 20 ) 


There are several ways of obtaining the right-hand side of (19). One could 
perform a quadrature for each k. We chose to sample 0(r) at N+l 
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points, approximate 0(r) by an N + 1 point approximation, and explicitly 
evaluate the required inner products. Equations (19) and (20) form an N x 
N system, which is tridiagonal plus a full last row. Special solvers may be 
used to solve such a system. 

Example 2 : 

Consider Bessel's equation of order 0. 

Lf(r) ■ 0 , L = [d(rd/dr)/dr + r] (21a) 

subject to either of: 

f(0) - J (0), or f(l) - J (1) (21b) 

o o 

We seek a (N+2) point solution in the form, 

N+l 

V r) < 2 » 

j“0 

A A 

(N + 2) equations for the (N + 2) unknowns f ... f^ + ^ are required. N 

equations are obtained from 

(Lf N+2 ,\ (r) r(l-r) 2 ) - 0 , k - 0, ... , (N - ) (23) 

To satisfy the boundary conditions, we set 

N+l 

f N+2 (0) " E f i V°> “ J o (0) (24a) 

j*0 J J 

N+l 

f N+2 (1) = E f i H 1 ( 1 > " V 1 ) (24b) 

j =0 J J 

A 

Thus, we obtain a matrix system for _f . This matrix is sparse, having 

entries along nine of its diagonals, as well as its last two rows. Again, we 

can utilize special matrix solvers. 
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III. Turbulent Flame Propagation 

A. Physics of Turbulent Flames 

The propagation of a flame front into a premixed fuel-oxidant mixture is 
important in a number of situations, including some types of gas turbine and 
internal combustion engines. The stability of such flames and the issues of 
ignition and extinction are other important issues in turbulent flame propa- 
gation. The importance of turbulent flame propagation and its difference from 
the axisymmetric jet problem described above made this problem our choice as 
the second case to be simulated. 

In all cases of interest in applications, the flame front is thin with 
respect to most of the other length scales of importance in the problem. In 
particular, actual flames are always thin with respect to the largest length 
scales of the turbulence. Whether they are smaller than the finest scales of 
the turbulence or not depends on such parameters as the Damkohler and Reynolds 
numbers. Both cases are of interest and importance and have been considered. 
The approaches to simulating the two cases differ somewhat. 

If the flame is thick with respect to the smallest scale of the turbu- 
lence (the Kolmogoroff scale), the appropriate approach is full simulation in 
which the chemical equations are solved together with the Navier-Stokes equa- 
tions; all of the details are then resolved on a single grid. This is very 
similar to the approach that has been used in nonreacting flows. However, 
with the computer resources available today, the parameter range accessible to 
this method is rather limited. 

On the other hand, if the flame is thin with respect to the Kolmogoroff 
scale, it is more efficient to use an approach in which the outer flow and the 
flame are computed on different scales and the solutions are matched. In this 
approach, the flame appears to be a thin discontinuity to the outer flow. The 
flame is computed on a smaller scale in which the flow appears laminar but may 
be unsteady, stretched, and/or curved. In this method, effort can be saved by 
precomputing the properties of the flame. This approach requires development 
of a method of dealing with discontinuities, as previous simulations of 
nonreacting turbulent flows have not required this feature. 

Because it offers the possibility of handling high Damkohler numbers, we 
expect to use the second approach but are reserving the possibility that the 
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first approach will be needed in the future. The method will be used to study 
the geometry of the flame front and its average speed of propagation as a 
function of various parameters (including the laminar flame speed and the 
length scales and intensity of the turbulence). 

B. Laminar Flame Calculations 

We have investigated several algorithms for solving the laminar flame 
problem in order to gain experience solving the equations, including solving 
both the unsteady and steady problems. We used a global, single-step kinetics 
model and assumed Lewis number unity. A test problem from Peters and Warnatz 
[3] served as the test. The nondimensional equations listed below were used 
to describe the laminar flame. 


l£ + ifiH = 0 

3t 3x 


where W 


3£U . 

, 2 
3pu 

3P (1) 

3t 

3x 

3x 

3pT 

3puT 

2 

3 T 

-Z— + 
3t 

3x 

- d + W 

2 

3pY 

3puY 

3x 

2 

3 Y 

+ 

3t 

3x 

= D — r ” W 
2 

P (°) 

- pT 

3x 

= constant 

T)/(l-ct 

C 1 — T ) ) ) 

is the reaction 


6 is the nondimensional activation energy, 
B is the Damkohler number, and 
a is the nondimensional heat release# 


For the unsteady problem, we used numerical techniques which could be 
extended to the multidimensional turbulent problem. Specifically, we used 
second-order central finite differences for the spatial derivatives and 
second-order Adams-Bashforth for the time integration# The unsteady problem 
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was solved with and without the convection terms in the energy equation. With 
the convection terms included, the problem is fairly difficult. This is 
discussed below. 


Solving the equations without the convection terras was fairly easy. In 
the standard method, the troublesome terms are removed by an appropriate 
coordinate transformation: 

x 

t = t, x . f pdx 

o 

where the overbar signifies the transformed coordinates. Unfortunately, this 
transformation cannot be used in higher dimensions. For the constant pressure 
case, the resulting equations are greatly simplified: 


lei 

at" 


th 



+ w 


M. 

at 



- w 


Using these equations, the internal structure of the flame is fairly easy to 
calculate, but the flame speed is extremely sensitive to the acccuracy of the 
calculation. The flame structure and speed from sample calculations are shown 
in Figs. 1—4. The structure appears well resolved with either 200 or 400 
points. However, the flame speed oscillates and the amplitude of oscillation 
depends on the number of points. The oscillations are due to small changes in 
the chemical source term as it propagates over the grid. Since the flame 

speed is dependent on the total heat release in the flame, these small changes 
are significant enough to cause the oscillations. The sensitivity of the lam- 
inar flame calculation to the source term indicates that special techniques 
such as grid refinement or grid adaptation will be very helpful. 

We also worked on solution techniques for the steady equations. The 

approach was to recast the constant-pressure problem as a system of ODEs and 

use standard subroutines to solve the system. The classical cold boundary 

problem now appears because the steady problem is actually ill— posed. 

First, we used a shooting technique. This involved guessing the initial 

temperature gradient and the flame speed, and integrating the equations 

through the flame. Correct initial guesses result in zero temperature gra- 
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dient at the end of the flame. As Is usual for stiff problems, shooting did 
not work very well. The problem could be written: 


Solve: 


Initial Conditions: 


g l = g 2 
§2 = P U S L ~ w 


gj(0) = 0 


g 2 (°) 


small initial guess 


where g^ = T 

g 2 = dT/dx 

W = B(l-T) exp(-f}( 1-T) / (l-a( 1-T) ) ) 

We had more success by posing the problem as a boundary-value problem and 
using a subroutine written by V. Pereyra called DVCPR. In this case, the 
flame-speed eigenvalue was treated as an unknown and another equation was 
added to the system. 

Define: Solve: Boundary Conditions: 



Adding the extra equation not only gives the desired solution without 
iteration, but it also allows an extra boundary condition to be enforced. 
This is highly desirable. This method worked very well, and the results were 
consistent with the Peters test problem results. For example, with a Damkoh- 
ler number of 50 and a nondimensional activation energy of 10, the flame 
velocity was 0.917, while the values reported by others ranged from 0.894 to 
0.918. 


13 


C. Evaluation of Numerical Methods 


We are restricting attention to low-speed flows to allow computations to 
be made without modeling. This is a common approach (Merkle and Choi [5]), so 
we shall list the nondimensionalized equations without development: 



|fi- + 

3t 

3pu 
3t + 

3 P u i u i 

3pT 

8p V 

3t 

3x j 

9pY ^ 

8p V 

3t 

3x 

j 


3pUj 

3x, 


3P 


( 1 ) 


3x, 


2 

3 T 


th 3x^ 3x^ 


2 

3 Y 

M 3x 3x 

j j 




+ W 


- W 


(25) 


P 


( 0 ) 


= pT 


constant 


where : 



°M 

W 


shear stress, 
thermal diffusivity, 
mass diffusivity, and 
chemical source term, f(Y,T). 


To find an effective numerical solution technique for these equations, we face 
two problems. Both arise from the decoupling of the thermodynamic pressure, 
P^®\ from the kinematic pressure, P^\ in the low Mach number approxima- 
tion and the density variation. 


One problem is best examined by combining the energy and the state equa- 
tions. The result is 


3u 

j 


—— D + W 

p (0) th 3x j 3x j 


(26) 


which clearly shows the main effect of the combustion on the velocity field. 
The problem arises from the lack of a time derivative in this equation. We 
have approached this problem in two ways. 
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In the first approach, we use the idea from incompressible flow that this 
equation is enforced via a Poisson equation for the pressure. Thus we take 
the divergence of the momentum equation and attempt to enforce Eq. (26) in the 
resulting equation. This approach could take several forms. We are examining 
the idea of decomposing the velocity field into two parts: one from the energy 
equation and one that is divergence-free. We are also looking at using a 
global projection operation similar to a "weak" formulation of the momentum 
equation [7]. Work is just beginning on this approach, and we have not yet 
found a suitable method. 


In the second approach, we solve one equation in nonconservative form; 
either the energ> equation or the continuity equation could be treated this 
way. McMurtry et al. [4] solve the continuity equation in nonconservative 
form and use the energy equation to replace the velocity divergence. Numer- 
ically solving an equation in nonconservative form requires extra care to 
assure the global conservation of certain quantities such as energy. In the 
low Mach number equations, we have the added difficulty that kinetic and 
thermal /chemical energy must be independently conserved. We would prefer to 
avoid nonconservative schemes if a suitable scheme can be found. 


The second problem which arises in solving the low Mach number equations 
concerns the Poisson equation for the pressure. When we take the divergence 
of the numerical momentum equation, the following term appears: 

n+1 

(where 6 is a numerical operator) 


6pu 4 


fix. 


For incompressible flows, this term is zero by continuity, 
continuity for this term in the low Mach number equations, 
approximate the term: 


If we use 
we have to 


- i£ 

fit 


n+1 


McMurtry uses this approach and approximates this terra by backward differ- 
encing. While this approach works, it would be better to avoid a numerical 
approximation of a time derivative in the Poisson equation. 

Finally, we note that these problems can be avoided completely by not 
using the low Mach number approximation. This would result in a new set of 
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problems. For example, in low-speed flows the time-step limitation can be 
severe. One approach used in low-speed compressible flows is to use acoustic 
subcycling and filtering. This approach is used by O'Rourke [6] and involves 
resolving the sound waves by integrating a subset of the equations over many 
small time steps during one regular time step. While this may work, it is too 
inefficient for our requirements. 

D. Two-Dimensional Code Development 

We are now developing a code for 2-D turbulent reacting flows. It is 
based on McMurtry's method, but this will be changed when we develop something 
better. The code is for premixed flames and uses single-step global kinetics. 
It is a full simulation code and thus has no modeling. In the streamwise 
direction we are using finite differences, because there are inflow and out- 
flow boundary conditions. In the spanwise direction, the boundaries are 
periodic, so a Fourier spectral method is used. The variable density Navier- 
Stokes solver and Poisson pressure solver are completed and in the testing 
phase. Solvers for the energy and species equations remain to be written but 
should be straightforward. 
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